function [lambda,weight] = quadpts3face(order)
%% QUADPTS3FACE quadrature points on a tetrahedron's faces
%
% [lambda,weight] = quadpts3face(order) return quadrature points with given
% order (up to 4) in the barycentric coordinates for faces
%
% The output lambda is a matrix of size nQ by 3, where nQ is the number of
% quadrature points. lambda(i,:) are the four barycentric coordinates of the
% i-th quadrature point and lambda(:,j) is the j-th barycentric coordinate
% of all quadrature points on each face. 
% The x-y coordinate of the p-th quadrature point
% can be computed as 
%
%     pxyz = lambda(p,1)*node(elem(:,1),:) ...
%         + lambda(p,2)*node(elem(:,2),:) ... 
%         + lambda(p,3)*node(elem(:,3),:) ... 
%         + lambda(p,4)*node(elem(:,4),:);
%
% The weight of p-th quadrature point is given by weight(p). See
% verifyquadpts for the usage of qudrature rules to compute integrals over
% triangles.
% 
% For the p-th quadrature point, it is on f-th face of the tetrahedron
% nQuadPerFace = nQ/4; f = ceil(p/nQuadPerFace);
% A sample of how to use this compute 
% \sum_{K}<phi_i, phi_j>_{\partial K}
%
% [lambdaAllFace, wAllFace] = quadpts3face(2);
% nQuadAllFace = size(lambdaAllFace,1);
% nQuadPerFace = nQuadAllFace/4;
% S = sparse(nDoF, nDoF);
% for i = 1:4
%     for j = 1:4
%         Sij = 0;
%         for p = 1:nQuadAllFace
%             Sij = lambdaAllFace(p, i)*lambdaAllFace(p, j);
%             f = ceil(p/nQuadPerFace); % current face
%             Sij = Sij + wAllFace(p)*Sij*faceArea2elem(:,f);
%         end
%         ii = elem2dof(:,i);
%         jj = elem2dof(:,j);
%         S = S + sparse(ii,jj,Sij,nDoF,nDoF);
%     end
% end
%
% References: 
%
% * Jinyun, Y. Symmetric Gaussian quadrature formulae for tetrahedronal
% regions. Comput. Methods Appl. Mech. Engrg.. 43(3):349--353, 1984.
%
% See also quadpts, quadpts3
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details. 

if order>4
    order = 4;
end
switch order
    %% Order 1, nQuad 1
    case 1 
        lambda = [0, 1/3, 1/3, 1/3;... %face 1
            1/3, 0, 1/3, 1/3; ...      %face 2
            1/3, 1/3, 0, 1/3; ...      %face 3
            1/3, 1/3, 1/3, 0];         %face 4
        weight = [1,1,1,1];
    %% Order 2, nQuad 3
    case 2
        lambda = [0, 2/3, 1/6, 1/6; ...%face 1
            0, 1/6, 2/3, 1/6; ...%face 1
            0, 1/6, 1/6, 2/3; ...%face 1
            2/3, 0, 1/6, 1/6; ...%face 2
            1/6, 0, 2/3, 1/6; ...%face 2
            1/6, 0, 1/6, 2/3; ...%face 2
            2/3, 1/6, 0, 1/6; ...%face 3
            1/6, 2/3, 0, 1/6; ...%face 3
            1/6, 1/6, 0, 2/3; ...%face 3
            2/3, 1/6, 1/6, 0; ...%face 4
            1/6, 2/3, 1/6, 0; ...%face 4
            1/6, 1/6, 2/3, 0];%face 4
        weight = [1/3, 1/3, 1/3, ...%face 1
            1/3, 1/3, 1/3,...%face 2
            1/3, 1/3, 1/3,...%face 3
            1/3, 1/3, 1/3];%face 4
    %% Order 3, nQuad 4
    case 3
        lambda = [0, 1/3, 1/3, 1/3; ...
            0, 0.6, 0.2, 0.2; ...
            0, 0.2, 0.6, 0.2; ...
            0, 0.2, 0.2, 0.6; ...
            1/3, 0, 1/3, 1/3; ...
            0.6, 0, 0.2, 0.2; ...
            0.2, 0, 0.6, 0.2; ...
            0.2, 0, 0.2, 0.6; 
            1/3, 1/3, 0, 1/3; ...
            0.6, 0.2, 0, 0.2; ...
            0.2, 0.6, 0, 0.2; ...
            0.2, 0.2, 0, 0.6;
            1/3, 1/3, 1/3, 0; ...
            0.6, 0.2, 0.2, 0; ...
            0.2, 0.6, 0.2, 0; ...
            0.2, 0.2, 0.6, 0];
        weight = [-27/48, 25/48, 25/48, 25/48, ...
            -27/48, 25/48, 25/48, 25/48,...
            -27/48, 25/48, 25/48, 25/48,...
            -27/48, 25/48, 25/48, 25/48];
    %% %% Order 4, nQuad 6
    case 4
        lambda = [0, 0.108103018168070, 0.445948490915965, 0.445948490915965; ...
            0, 0.445948490915965, 0.108103018168070, 0.445948490915965; ...
            0, 0.445948490915965, 0.445948490915965, 0.108103018168070; ...
            0, 0.816847572980459, 0.091576213509771, 0.091576213509771; ...
            0, 0.091576213509771, 0.816847572980459, 0.091576213509771; ...
            0, 0.091576213509771, 0.091576213509771, 0.816847572980459; ...%end of face 1
            0.108103018168070, 0, 0.445948490915965, 0.445948490915965; ...
            0.445948490915965, 0, 0.108103018168070, 0.445948490915965; ...
            0.445948490915965, 0, 0.445948490915965, 0.108103018168070; ...
            0.816847572980459, 0, 0.091576213509771, 0.091576213509771; ...
            0.091576213509771, 0, 0.816847572980459, 0.091576213509771; ...
            0.091576213509771, 0, 0.091576213509771, 0.816847572980459; ...%end of face 2
            0.108103018168070, 0.445948490915965, 0, 0.445948490915965; ...
            0.445948490915965, 0.108103018168070, 0, 0.445948490915965; ...
            0.445948490915965, 0.445948490915965, 0, 0.108103018168070; ...
            0.816847572980459, 0.091576213509771, 0, 0.091576213509771; ...
            0.091576213509771, 0.816847572980459, 0, 0.091576213509771; ...
            0.091576213509771, 0.091576213509771, 0, 0.816847572980459; ...%end of face 3
            0.108103018168070, 0.445948490915965, 0.445948490915965, 0; ...
            0.445948490915965, 0.108103018168070, 0.445948490915965, 0; ...
            0.445948490915965, 0.445948490915965, 0.108103018168070, 0; ...
            0.816847572980459, 0.091576213509771, 0.091576213509771, 0; ...
            0.091576213509771, 0.816847572980459, 0.091576213509771, 0; ...
            0.091576213509771, 0.091576213509771, 0.816847572980459, 0]; %end of face 4
        weight = [0.223381589678011, 0.223381589678011, 0.223381589678011, ...
            0.109951743655322, 0.109951743655322, 0.109951743655322, ...
            0.223381589678011, 0.223381589678011, 0.223381589678011, ...
            0.109951743655322, 0.109951743655322, 0.109951743655322, ...
            0.223381589678011, 0.223381589678011, 0.223381589678011, ...
            0.109951743655322, 0.109951743655322, 0.109951743655322, ...
            0.223381589678011, 0.223381589678011, 0.223381589678011, ...
            0.109951743655322, 0.109951743655322, 0.109951743655322];
    
end
